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ABSTRACT 

This paper presents a technique for quantifying the wear or damage of gear teeth in a transmission system. 
The procedure developed in this study can be applied as a part of either an onboard machine health- 
monitoring system or a health diagnostic system used during regular maintenance. As the developed 
methodology is based on analysis of gearbox vibration under normal operating conditions, no shutdown or 
special modification of operating parameters is required during the diagnostic process. 

The process of quantifying the wear or damage of gear teeth requires a set of measured vibration data and a 
model of the gear mesh dynamics. An optimization problem is formulated to determine the profile of a 
time-varying mesh stiffness parameter for which the model output approximates the measured data. The 
resulting stiffness profile is then related to the level of gear tooth wear or damage. 

The procedure was applied to a data set generated artificially and to another obtained experimentally from a 
spiral bevel gear test rig. The results demonstrate the utility of the procedure as part of an overall health- 
monitoring system. 

Keywords: Health monitoring; Gears; Maintenance; Wigner-Ville distribution; Time-frequency analysis; 
Optimal tracking 


1. INTRODUCTION 

In the last two decades, with demands for higher operating speeds and greater load capacity, premature 
failures in high-performance turbomachinery have often resulted in enormous financial losses and, at times, 
catastrophic consequences. In aeronautical applications, where both weight and efficiency are pushed to their 
design limits, the prevention and management of premature equipment failures is a vital part of the 
maintenance program. Current onboard condition-monitoring systems for gas turbine engines often fail to 
provide sufficient time between warning and failure for safety procedures to be implemented. On the other 
hand, inaccurate interpretation of operating conditions may result in false alarms and unnecessary repairs and 
downtime. The early detection of incipient failure in a mechanical system is of great practical importance as 
it permits scheduled inspections without costly shutdowns and indicates the urgency and locations for repair 
before a system incurs catastrophic failure. 

Some success has been achieved in identifying damage in a gear transmission system by using a joint 
time-frequency analysis known as the Wigner-Ville distribution (WVD) technique (Boashash and Black, 
1987; Choy et al., 1994a,b; Claasen and Mecklenbrauker, 1980). The approach is to use statistical pattern 
recognition to match the WVD signature patterns of damaged gears with standard patterns stored in a data 
base. Although the WVD technique is useful for determining the type and location of the damage, it is not 
much help in quantifying the level of damage. Damage quantification would logically be the next step in 
failure prediction; however, no explicit attempts at damage quantification have previously appeared in the 

literature. . 

This paper presents a new technique for processing vibration data to quantify the level of damage in a gear 

transmission system. The technique consists of a nonlinear numerical optimization in the form of an 
“optimal tracking” problem (Sage, 1968; Lewis, 1986). The optimization uses a dynamic model of the 
gear mesh and forms an estimate of the time- varying mesh stiffness that best corresponds to the given set of 
vibration data. The utility of the technique relies on the relationship between the wear or damage of a gear 
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tooth and the change in stiffness of the mesh during a given tooth pass cycle. An analysis of this 
relationship demonstrates that the perturbation of the stiffness function from the nominal profile can be 
used to quantify the level of damage. 

The optimal tracking technique for estimating the perturbation of the mesh stiffness was tested in two 
settings. First, it was tested on a set of fictitious data generated by computer simulation of a one-degree-of- 
ffeedom mechanical system with time-varying stiffness. The solution of the optimal tracking problem 
matched very closely the actual stiffness profile used in the model generating the data. Then, the technique 
was tested on a set of experimental data from a gear test rig, but still assuming the one-degree-of-ffeedom 
model. Despite the simplicity of the model the stiffness profile obtained was shown to be useful in 
correlating to the level of damage of the gear transmission system. 

The paper is organized as follows: Section 2 presents the system model and formulates the optimal 
tracking problem. Section 3 outlines the numerical solution procedure for the nonlinear optimization. 
Section 4 presents and interprets the results of the optimization and discusses the next steps to be taken in 
developing a comprehensive failure-prediction procedure. 


2 . OPTIMAL TRACKING PROBLEM 
2*1 SYSTEM MODEL 

The system considered in this study consisted of a small pinion in mesh with a larger gear. A simple 
model of this system has the two gear masses connected by a spring and a damper. The larger gear is much 
heavier than the pinion; hence, it is assumed to be rigid, so that all relative motion between the two is 
attributed to the motion of the pinion. Then, the equation of motion of the pinion takes the form 

mx + cx + k(t)x = 0, (1) 

where m is the mass of the pinion and k(t) and c are the stiffness and damping of the mesh. The mesh 
stiffness is not constant but is nominally a periodic function of the gear angle, with each period 
corresponding to one tooth pass. The high points on the periodic stiffness function correspond to gear 
angles where two pairs of gear teeth are in contact, and the low points correspond to angles where only one 
pair is in contact. 

It has been found in experiments on gearbox vibrations (Choy et al., 1994b, 1995) that the gear mesh 
stiffness changes with the wear, pitting, or fracture of the gear teeth. Such changes in the gear mesh 
stiffness inevitably lead to changes in the vibration signatures of the mechanical system. The objective of 
the optimal tracking procedure developed in this study is to reconstruct the true stiffness profile for a 
damaged gear tooth from the measured vibration. That is, the objective is to determine the function k(t) that 
would result in the measured vibration according to the system model (1). 

The true stiffness profile can be expressed as the sum of a constant (time averaged) component, a nominal 
periodic component, and a perturbation resulting from gear wear or damage. Accordingly, the system model 
(1) is written as 


mx + cx + 


^ave aperiodic ^ ^perturb 


(/)]* = o. 


(2) 


or 

x + — x + £2 2 x = u(t)x, (3) 

m 

where Q 2 = k^^fm and u(t) represents the total time- varying component of the stiffness divided by the 
pinion mass. By defining the variables 


x { = i, x 2 = x, 


(4) 


2 



the system model can be written in the state- variable form 


x, = — —x. -£2 2 x 2 +u(t)x 2 
m 

(5) 

x 2 = x \ 


with the given initial conditions 


X 2 ^0 ^ * 0 ' 


( 6 ) 


2.2 OPTIMIZATION PROBLEM 

Suppose that a data set corresponding to the vibration of the pinion is collected over the interval 
Let the function describing the data set be denoted as x 2 (r), since it corresponds to the modeled variable 

x 2 (/). The objective is to determine a reasonable time-varying stiffness component u(t ) for which the model 
output *2 (0 approximates the measured data x 2 (t). 

A diagram depicting the functional objective is shown in Fig. 1(a). In the figure u(t ) is depicted as an 
input to be chosen so that the error e(t) will be small for all time. Note that this problem has the form of a 
tracking problem, where the control input of a system is designed so that the system output follows a 
prescribed reference function. Such a problem may be approached by using the standard techniques of 
optimal control theory (Sage, 1968; Lewis, 1986). In particular, the “design” of a suitable function u(t) 
may be achieved by minimizing the quadratic cost functional 

J = ^P f [x 2 (t f )-x 2 (t f )f +^{h[x 2 (t)-x 2 (t)f +p 2 u 2 (t)yt, (7) 

f 0 


where p x , P 2 , and /L are cost-function weighting coefficients. This form of the cost functional penalizes the 
energy in the error between the modeled output and the measured data. It also penalizes the use of too large 
a stiffness perturbation function in order to avoid singularity in the solution. 

In the optimal tracking problem the system dynamic equations (5) are treated as equality constraints 
imposed in the optimization of the cost (7). As such, they are appended to the cost function by using time- 
varying Lagrange multipliers Aj(r) and ^(r). These Lagrange multipliers are themselves governed by 
differential equations called the costate equations. The costate equations together with the state equations of 
the system model form a two-point boundary value problem (TPBVP) (Sage, 1968; Lewis, 1986). The 
TPBVP equations are 


(State equations) 


x l =- — x l -Q 2 x 2 + u(t)x 2 

m 


x 2 =x x 


( 8 ) 
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(Costate equations) 


(9) 


Aj — ~ A2 H Aj 

m 


X 2 =Q 2 X i -u(t)X i -p l [x 2 (t)-X 2 (t)] 


(Stationarity condition) 0 = AjX 2 + /J 2 u(r) (10) 


(Endpoint conditions) x x (/ Q ) = Xq , * 2 (*o ) = *0 (1 1) 

A,('/> = 0 ’ X 2 (t f ) = p f [x 2 (t f )-x 2 (t f )\. (12) 

The TPBVP (8-12) represents a set of necessary conditions for u(r) to be the solution of the optimal 
tracking problem. The TPBVP consists of a set of four coupled differential equations (8-9), together with 
an algebraic relation (10), and some endpoint conditions (1 1-12) at both t 0 and t r Notice that the TPBVP is 
nonlinear: the unknown function u{t) multiplies other independent variables in the differential equations. 


3. NUMERICAL SOLUTION PROCEDURE 

The nonlinear TPBVP (8-12) is solved by an iterative procedure. A complete and general derivation of the 
procedure is given in Sage (1968) and Dyer and McReynolds (1970). Some of the salient points are outlined 
below for convenience. 


3.1 SUCCESSIVE SWEEP METHOD 

Solving the nonlinear TPBVP requires an iterative method. Although several approaches are possible, a 
common and useful one is to begin with an initial guess u\t) and use it to integrate the nonlinear state 
equations (8) forward in time starting from the initial conditions (11) to determine the nominal state 
functions x®(t) and x 2 (f). Then, starting from the final conditions (12), integrate the nonlinear costate 

equations backward in time to determine the nominal costate functions Aj°(r) and A^( t ). The nominal 
functions u°(r), * 2 (r), A^(0, and A^(f) then satisfy all the TPBVP equations except the 

stationarity condition (10). 

The nominal state, costate, and input functions must be iteratively updated, so that they will eventually 
satisfy all the nonlinear TPBVP equations, including the stationarity condition. Each update is 
accomplished by solving a linearized version of the TPBVP. A standard method for doing this is known as 
the sweep method, whereby a linear relationship between the state and costate functions is assumed. Then, 
the linear TPBVP degenerates into a set of ordinary differential equations with endpoint conditions at the 
final time only. These are solved by a straightforward numerical integration. In the case of the optimal 
tracker these ordinary differential equations take the form of the coupled Riccati equations 
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(13) 


P\ \ = 2 — Pn ~Pn + Pi\ 


.2 *L 


02 


P \2 = — Pn ~ P 22 ~ P\ 1 

m 


-Q 2 +u(t)- — X\X 2 

0 2 


x\ 

+ PnPl2 ~n~ 

P 2 


P 22 - _ 2pi2 


-n 2 +u(t)-—x x x 2 

P 2 


r 2 

, 2 x 2 
+ p i2ir~ 
02 


f 


with endpoint conditions 


P\\(. t f)-P\2^ t f^>~^' P22^f^~Pf’ 


(14) 


together with the auxiliary linear equations 

_2\ 


= *1 


V 


x 2 

c+pn h 


-h 2 -p n e^[l ] x 2 +P 2 u(t)} 

P 2 


(15) 


h 2 = -hi 




--7r[*-l x 2 + fa u( J)]{P\2 x 2 + M ) 

P 2 


with the endpoint conditions 


h l (t f ) = h 2 (t f ) = 0. 


(16) 


Note that the x and A variables in the differential equations (13) and (15) represent the given nominal 
functions. (The zero superscripts have been omitted for convenience.) They are simply treated as time- 
varying coefficients in the numerical integration of the differential equations. The solutions of 
equations (13) and (15) are then used to compute the corrections to the nominal state, costate, and input 
functions. This computation requires yet another numerical integration, this time of the linearized state 
equations 


Aij = Axj 


x 2 

c ; p "hj 


-I- Ax. 


-Q 2 +u(t)-—k l x 2 -Pi 2 -ir 
p 2 Pl m 


P 2 P 2 


Ax 2 = Axj 


(17) 


with the zero initial conditions 


Ax,(f 0 ) = Ax 2 (/ 0 ) = 0. 


(18) 


Finally, the update of the nominal control is computed as 

_1 

02 


Au = ^[^x 2 +p 2 u(t)]-j-[?i l Ax 2 + x 2 (p n tei+Pi 2 *x 2 +hi)\ 
02 


(19) 
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where € is the step size, and the new nominal control is given by 


u M (t) = «'(/)+ A«' (/)• (20) 

(The superscripts i and i+1 denoting the iteration number have been reinserted in equation (20).) The 
procedure is repeated until the nominal functions converge to a solution. 

The real scalar e e [0,1] in equations (15), (17), and (19) is used as a “step size” parameter. Using a 
smaller value of e tends to decrease the magnitudes of the corrections, thereby improving the stability of the 
iterative procedure but slowing the convergence to the solution. Using a larger value of € has the opposite 
effects. 


3.2 NUMERICAL DETAILS 

The choices of the cost-function weighting coefficients ft, ft, and ft are important for effective numerical 
optimization. The parameter ft defines the penalty on the difference between the calculated and reference 
vibration signals. Since the goal is to minimize the difference between the calculated and tracked vibration 
signals, a large weighting coefficient ft should be chosen. The parameter ft defines the penalty on the 
function u(r). Generically speaking, ft should impose a lighter penalty on u(t) than ft imposes on the 
tracking error. Note also that the choice of units has an effect on the appropriate relative sizes of ft and ft. 
In the examples studied the numerical values of u(t) are considerably larger in magnitude than those of a 
reasonable vibration-signal error; therefore, even if equal weighting between error and control were desired, 
ft should be chosen to be considerably smaller than ft. An inappropriately large choice of the parameter ft 
would make the cost function almost unchanged from one iteration to the next. Thus, a small constant 
value was chosen for the parameter ft. The parameter 1 3 f defines the penalty for the error at the final time. 
If Pf is too small, a large vibration error at the final time will result. 

By following these general guidelines the optimization algorithm described in the previous section was 
realized in a computer program. The equations were integrated with a seventh-order Runge-Kutta-Fehlberg 
method. A summary of the programming steps is given below (fig. 1(b)): 

0. Set i = 0 and take the initial guess «®(r) for the function u(t) to be zero. 

1 . Using the function u l (t) from the previous step, integrate the state equations (8) forward in time. 

Calculate the resulting cost function J l . 

2. Integrate the costate equations (9) backward in time. 

3. Use the computed state and costate variables as time-varying coefficients in the integration of the 

Riccati equations (13) along with the auxiliary equations (15) backward in time. 

4. Integrate the linearized state equations (17) forward in time. Using the linearized stationarity condition 

(19), calculate the correction A u*(t) to the nominal function u l (t) and hence the updated function 

a l+ *(r). Also, calculate the new cost function / l+ *. 

5. Make decisions about the continuation of the optimization procedure and the choice of the parameters: 

a. If the difference between the calculated and tracked vibration signals is small, the optimization 
procedure is finished. 

b. If the difference / l+ * - J l < 0 is large enough, repeat from step 1. 

c. If the difference 7* + * - J l < 0 is too small, increase the weighting ft and repeat from step 1 . 

d. If / 1+ * > A repeat from step 1 using a smaller value of the step size e. If this is not successful, 

increase the error weighting ft and repeat from step 1. 

Some comments should be made on step 5 of the numerical procedure. It was observed that for given 
values of weighting coefficients and the step-size parameter the optimization procedure converges to some 

value of the cost function. In this case the difference between the values of the cost functions - J 1 
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becomes negligible after some iterations. This means that the cost associated with the control w(f) is 
becoming dominant. Therefore, it makes sense to start a new iteration with an increased weight ft (i.e., 

imposing a higher penalty on the vibration error). 


4. DISCUSSION OF RESULTS 

To demonstrate the optimal tracking procedure described above, two numerical cases were used in this 
study. The first case was a numerical experiment in which the tracker was applied to a set of vibration 
signals generated numerically, assuming a given gear mesh stiffness profile. The mesh stiffness profile 
evaluated by the optimal tracker was compared with the original stiffness used in generating the vibration 
signal. Figure 2(a) shows the comparison between the vibration signal generated by a sinusoidal stiffness 
and that simulated by the optimal tracker. As shown in the figure the two vibration signals were very 
similar; the small difference between the two signals is given in Fig. 2(b). Figure 3(a) depicts the original 
gear mesh stiffness used and the stiffness evaluated by using the optimal tracker; the difference between the 
two stiffnesses is given in Fig. 3(b). The excellent agreement between the two stiffnesses in this numerical 
experiment has confirmed the applicability of the optimal tracking procedure in evaluating system stiffness 
changes from system vibration signals. However, this close resemblance between the generated and 
simulated signals was partly due to the original time signals being smooth, continuous, and harmonic 
without any substantial change in magnitude and phase over the gear revolution. To demonstrate the 
generality and limitation of the developed procedure, a set of experimental data taken from a test rig was 
used in the next case. 

The second case was based on the experimental data obtained from the spiral bevel gear test ng shown in 
Fig. 4. The primary purpose of this rig is to study the effects of gear tooth design, gear materials, and 
lubrication types on the fatigue strength of aircraft-quality gears (Zakrajsek et al., 1994). Because spiral 
bevel gears are used extensively in helicopter transmissions to transfer power between nonparallel 
intersecting shafts, using this fatigue rig for diagnostic studies is extremely practical. Vibration data from 
an accelerometer mounted on the pinion shaft bearing housing were captured by using a personal computer 
with an analog-to-digital conversion board and an anti-aliasing filter. The 12-tooth test pinion and the 36- 
tooth gear have the following characteristics: 0.5141 in pitch, 35° spiral angle, 1-in. face width, 90 shaft 
angle, and 22.5° pressure angle. The pinion transmits 720 hp at a nominal speed of 14 400 rpm. The test 
rig was started and stopped several times for gear damage inspection. The test was ended at 17.72 
operational hours when a broken portion of a tooth was found visually during one of the shutdowns. 

Figure 5(a) depicts the gear tooth after 6.5 hr of operation. Note that there is heavy surface pitting on one 
gear tooth with minor pitting on the next tooth. Figure 5(b) shows the time domain averaging, the 
frequency spectrum, and the joint time-frequency analysis using the Wigner-Ville distribution (WVD) 
(Boashash and Black, 1987; Choy et al„ 1994a,b, 1995; Claasen and Mechlenbrauker, 1980) of the 
accelerometer signal at 6.5 hr (Choy et al„ 1994a). Note that in Fig. 5(b) the time signal indicates a large 
vibratory signal during the engagement of the sixth and seventh teeth (damaged teeth), but the frequency 
spectrum, because of its averaging characteristics, shows very little change from the original signal (Choy 
et al., 1994a). The WVD begins to show a pattern of shifting of the major frequency component (at a mesh 
frequency of 2880 Hz) around the meshing of the sixth and seventh teeth. The WVD pattern in this case is 
very similar to those resulting from a short-term amplitude and phase change of a vibration signal (Choy 
et al., 1994a). Although it has been established by the authors in some previous publications (Choy 
et al., 1994a,b, 1995) that such damage in the gear can be identified by the WVD pattern recognition 
process, the level of the damage has not been addressed. A recent study by the authors has shown that wear 
and surface pitting of the gear tooth usually will result in a phase shift in the stiffness profile, without any 
significant change in the stiffness magnitude. Figure 6 shows the stiffness change in a gear mesh evaluated 
(Boyd and Pike, 1985) from gear tooth surface profile variations. Note in Fig. 6(b) that increasing surface 
profile variation increases the phase shift of the gear stiffness without changing the magnitude of the 
stiffness. 

Incorporating this constant gear mesh stiffness as an additional constraint, the optimal tracking procedure 
was applied to the experimental vibration signal (obtained from the bevel gear test rig at 6.5 hr as shown in 
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Fig. 5) to evaluate the corresponding gear mesh stiffness. To better evaluate the gear mesh stiffness, the 
time signal was filtered at a mesh frequency of 2880 Hz. Figure 7(a) shows the comparison between the 
unfiltered experimental signal and the optimal tracker simulation, and Fig. 7(b) shows the comparison 
between the filtered experimental signal and the tracker-simulated signal. Note that because of the 
substantial change of magnitude and phase of the time signal during the data acquisition period (one 
revolution of the gear), the accuracy in the simulated vibration is not as good as that in the numerical 
experiment (Fig. 2(a)). Figure 8 depicts the gear mesh stiffness evaluated by using the optimal tracker. 
Note that in the evaluated stiffness considerable phase shifts at several gear teeth resulted in the large 
variation in magnitude and phase of the vibration signal. At the location where pitting occurred (teeth 6 and 
7) the phase shift of the stiffness was more pronounced. By using the results from the evaluated mesh 
stiffness and the correlation of stiffness change with gear wear shown in Fig. 6(b), the gear damage can be 
estimated. 

5. CONCLUSIONS 

This paper presents a unified approach to identifying and quantifying damage in a gear transmission 
system. The conclusions from this study are as follows: 

1 . The application of the joint time-frequency technique called the Wigner-Ville distribution provides the 
ability to identify the types and locations of the gear damage. 

2. The optimal tracker developed in this paper provides a very reasonable estimate of the stiffness change 
due to damage, which can be related to the level of gear damage. 

3. For vibratory signals with large changes in magnitude and phase angle the accuracy of the simulated 
signal from the optimal tracker may decrease. 

4. For a more accurate evaluation of system mesh stiffness an optimal tracker for the complete dynamic 
model of the system is needed. 
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Figure 1 . — Block diagram of tracking problem. 
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Figure 5. — Damaged gear tooth and corresponding analysis, (a) Damaged gear 
tooth at 6.5 hr. (b) Time signal, WVD, and frequency spectrum of vibration 
signal at 6.5 hr. 
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Figure 6. — Simulated gear tooth wear and corresponding gear stiffness, (a) Simulated 
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calculated simulated tooth wear. 
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Figure 8. — Tracker-evaluted mesh stiffness for spiral bevel gear at 6.5 hr. 
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of the procedure as part of an overall health-monitoring system. 
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